Wind and boundary layers in Rayleigh-Benard convection. 
II. Boundary layer character and scaling 
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The scaling of the kinematic boundary layer thickness A„ and the friction factor C/ at the top- 
and bottom- wall of Rayleigh-Benard convection is studied by Direct Numerical Simulation (DNS). 
By a detailed analysis of the friction factor, a new parameterisation for C/ and A„ is proposed. The 
simulations were made of an L/H — 4 aspect-ratio domain with periodic lateral boundary conditions 
at Ra = {10^, 10^, 10^, 10*} and Pr = 1. A continuous spectrum, as well as significant forcing due 
to Reynolds stresses clearly indicates a turbulent character of the boundary layer, while viscous 
effects cannot be neglected judging from the scaling of classical integral boundary layer parameters 
with Reynolds number. Using a conceptual wind model, we find that the friction factor C/ should 
scale proportional to the thermal boundary layer thickness as C/ oc Ae, while the kinetic boundary 
layer thickness Au scales inversely proportional to the thermal boundary layer thickness and wind 
Reynolds number Xu oc Ag^Re^^. The predicted trends for C/ and \u are in agreement with DNS 
results. 

PACS numbers: 44.25.-|-f, 47.27.ek, 47.27.eb, 47.27.te 

Keywords: Rayleigh-Benard convection, wind, DNS, boundary layer, friction factor, turbulence 



I. INTRODUCTION 



The structure of the boundary layer is of great impor- 
tance for understanding the turbulent heat transfer char- 
acteristics of Rayleigh-Benard convection. Inherently un- 
stable due to buoyancy effects, the thermal boundary 
layer with thickness Ae is in a dynamic equilibrium of 
heating (cooling) by thermal diffusion and the detrain- 
ment (entrainment) of heat due to impinging and eject- 
ing thermals at the bottom (top) plate. This process 
creates large temperature gradients across the boundary 
layer, thereby enhancing the heat transfer through the 
wall and thus the Nusselt number Nu. Next to a thermal 
boundary layer, one can identify a kinematic boundary 
layer with thickness A„, associated with the velocity field. 
Depending on the Prandtl number Pr = i/n^^, which is 
the ratio between the kinematic viscosity i/ and thermal 
diffusivity k, the kinematic boundary layer can be nested 
inside the thermal boundary layer or vice versa, which in- 
fluences the effectiveness of the heat transfer as a function 
of the Raylcigh number Ra. The Raylcigh number Ra is 
defined as Ra = (3gAQH^{i'K)~^, where /? is the thermal 
expansion coefficient, g the gravitational constant, AG 
the temperature difference between the top and bottom 
plate and H the domain height. The scaling of Ae and 
A„ as a function of Ra and Pr are therefore of importance 
for proper prediction of the heat transfer. 

In the theory of Grossmann and Lohse the wind 



velocity U and the boundary layer thicknesses A^ and Ae 
are central parameters, which are used to estimate the 
dissipation rates of kinematic energy and temperature 
variance in the bulk and the boundary layers. In the 
theory, Ae and Au are defined as 



Ae oc i?/(2Nu), 
A„ cx HRc-^^^. 



(1) 
(2) 
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While ID) holds excellently, the correspondence of ^ 
with experiments [3, Q and simulations [1, Q is less 
satisfactory. Relation ([2]) can be obtained by non- 
dimensionalising the steady laminar two-dimensional 
Prandtl boundary layer equations @, 0|i by which ^ 
follows immediately. However, the measured Re depen- 
dence of A„ is much weaker than predicted by ^ (see 
also Fig. [!])■ It has been suggested that the difference 
is due to geometry effects (plate-filling vs. laterally 
restricted flow). 

In this paper, we argue that the disparity in the ex- 
pected and the observed scaling is because the top- and 
bottom boundary layers are not laminar, i.e. forcing due 
to Reynolds stresses cannot be neglected in the kinematic 
boundary layer. Consequently, the arguments leading 
to (HD do not hold. With a detailed DNS study of the 
momentum- and heat-budgets and the friction factor, 
and using the wind model of the accompanying paper 
0, we derive new parameterisations for A-^ and C/. 

A related question is whether or not the boundary lay- 
ers can be regarded as turbulent. The Reynolds number 
Re is too low to sustain a "classical" turbulent bound- 
ary layer (Re ~ 1500 at Ra ~ 10®), i.e. a boundary 
layer where the turbulence production due to shear is 
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in local equilibrium with dissipation. Hence, the general 
view is that the boundary layers are laminar, but time- 
dependent. Although time-dependence due to plume im- 
pingement and detachment prevents laminarity in the 
strict sense, the assumption could be justified if the 
plumes are passive with respect to the scaling of integral 
boundary layer parameters such as the friction factor C / 
and the kinematic boundary layer thickness Xu- Several 
other studies show that the friction factor scales simi- 
lar to a Blasius boundary layer [l^, [HI, However, 
the scaling of A„ does not comply with classical laminar 
scaling as discussed before. Furthermore, a recent 
study of time-spectra in the bottom kinematic bound- 
ary layer revealed that the spectra in the boundary layer 
and in the bulk were practically indistinguishable [Tsj , a 
strong indication of turbulence. In order to understand 
this dual behavior, we study several turbulence indica- 
tors for the boundary layers, such as the spectra, shape- 
and friction factor. 

The paper is outlined as follows. A brief summary of 
the code for direct simulation and symmetry-accounted 
ensemble-averaging is given in section [TTl The scaling 
of the boundary layer thickness, the velocity profile, the 
friction factor and the shape factor are studied in sec- 
tions nil A[ IIIIBI and IIIICI respectively. Then, we study 
the space and time spectra (section IIII D[) . In section 
nil El the mean momentum and temperature budgets in 
the boundary layers are studied to clarify the importance 
of fluctuations in the boundary layers. Using the results 
from the momentum budgets, the friction factor C/ is de- 
composed in a pressure and momentum-flux contribution 
in section UlI Fl This leads to the insight that the main 
contribution is by the pressure gradient. Using the con- 
ceptual wind model derived in the accompanying paper 
0, scaling laws for C/ and A„ are derived in section HVl 
As the results show that the flow has many typical fea- 
tures of turbulence but also of laminarity, the interpreta- 
tion of the results is discussed in section [V] Conclusions 
are drawn in section IVTl 



II. SIMULATIONS 

Direct simulation of Rayleigh-Benard convection has 
been performed at Ra = {10^ 10^, lO'^, 10*} and Pr = 1 
in a r = 4 aspect-ratio domain. The code is based on 
a second-order variance-preserving finite-difference dis- 
cretisation of the three-dimensional Navier-Stokes equa- 
tions and is fully parallellised. For all simulations, a grid 
was used sufficient resolution to resolve the smallest tur- 
bulent scales, i.e. the Kolmogorov scale rjK = {v^ /eY^'^ 
and Corrsin scale rjx = Pr Vk- The top and bot- 
tom wall are rigid (no-slip) and of fixed temperature. At 
the side domain boundaries, periodic boundary condi- 
tions are applied. For each Ra except the highest, 400 
independent realisations were obtained by performing 10 
independent simulations and sampling the velocity and 
temperature field roughly twice every convective turnover 



time. Because of a formidable computational require- 
ments for Ra =10*, we use this simulation only for the 
results of Fig. [1] and confine the wind-decomposed anal- 
ysis to the lower Ra cases, though without loss of gener- 
ality. 

Similar to domains confined by sidewalls, a wind 
structure develops also in domains with lateral periodic 
boundary conditions. However, here the wind structure 
can be located anywhere in the domain since it is not kept 
in place by sidewalls. To extract the wind, symmetry- 
accounted ensemble-averaging is used [l3 |. which aligns 
the wind structure in different realizations before averag- 
ing. In this way a wind structure can be identified unam- 
biguously for these domains, by which a decomposition 
in wind and fluctuations becomes possible. The result- 
ing average velocity and temperature (three-dimensional 
fields) are denoted respectively by Ui and Q. The tildes 
are used to distinguish the conditional average from the 
standard (long-time, ensemble or plane) average X which 
is a function of z only. The symmetry-accounted average 
can be interpreted exactly as classical Reynolds-averaged 
results. For further details we refer to the accompanying 
paper 

III. RESULTS 
A. Boundary layer thickness 

The thickness of the hydrodynamic and thermal 
boundary layers as a function of Ra is shown in Fig. 
[T] Here, All and Ae arc defined as the location of 
the maximum of mean squared horizontal velocity fluc- 
tuations u'u' and mean squared temperature fluctua- 
tions respectively. The approximate powerlaws are 
A„ = 0.5 Ra"° " and Ae = 2.33 Ra"°-^^ respectively, in 
good agreement with other simulations [1] and reason- 
able agreement with experiments [2[ (despite differences 
in aspect-ratio, geometry and boundary conditions). 

Also shown in Fig. [T] are the predictions of the bound- 
ary layer thickness ([1]), ^ from the Grossmann-Lohsc 
theory [l[, together with the DNS results. The thermal 
boundary layer thickness Ae is in good agreement with 
the simulations. The width of the kinematic boundary 
layer A„ does not agree so well with the G-L theory, as 
the actual exponent is —0.13 instead of —0.25 (where we 
have assumed free-fall scaling Re oc Ra^^^ for simplicity). 

Below we briefly recapitulate the arguments of [J 
leading to The starting point is the laminar two- 
dimensional Prandtl boundary layer equation [^, [l^ 

udxU + wdzU = vd^u. (3) 

1/2 

Upon substituting x — > Hx, z HUc ' z, u Uu 

I/O 

and w — > C/Re ' w, the equations become parameter 
independent as 

udxU + wdzU — d^u. (4) 
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FIG. 1: The thermal and kinematic boundary layer thickness 
Ae and Au as a function of Ra. The dashed and dash-dotted 
lines in the graphs are predictions by the G-L theory of Au 
and Ae, respectively. 



Neither this expression, nor the incompressibility condi- 
tions, nor the boundary conditions have an cxpUcit de- 
pendence on Re, so the solution has to be independent of 
Rc as well. Therefore, the flow pattern undergoes a sim- 
ilarity transformation, and the boundary layer thickness 
scales as \u/H cx Rc^^''^. This result is rigorous, pro- 
vided that ^ holds, i.e. that turbulent stresses do not 
play a role in the momentum budget. In section IIII El 
we show that forcing due to Reynolds-stresses cannot be 
neglected for the boundary layer equations so that the 
laminarity assumption does not hold. 



B. Velocity profiles 

The characteristic shape of the velocity profile can be 
obtained from the plane-averaged horizontal average ve- 
locity as u{z) ~ ^ (uu) A ■ Figure [2k shows these profiles 
for various Ra in plus units, i.e. scaled by the friction 
-^/t^^^Tp with ~ u/ut and ~ zur/v- 

as 



velocity u-r = \JtwI p with w+ ~ ujur and z+ 
Here, we define a typical wall shear-stress 



= pvdz \/ {uu) 



A\w 



(5) 



In Fig. the viscous sublayer relation u+ = z~^ is shown 
with a dashed line, and logarithmic scaling of the velocity 
profile results in a straight line. For a classical turbulent 
channel flow and constant-pressure boundary layer, the 
viscous region ends at y"*" = 5, the log-layer starts from 
~ 30 and the profiles will collapse onto a single uni- 
versal curve for all Re. Here the situation is completely 
different. First, in plus-coordinates the profiles do not 
collapse onto a single curve. Furthermore, the viscous 
region ends at approximately = 1, and the velocity 
reaches its maximum at z~^ « 10 at Ra = 10^. A re- 
gion where the velocity scales logarithmically is hard to 
distinguish, indicating the absence of an inner (constant 
stress) layer. 
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FIG. 2: The horizontal boundary layer velocity profile based 
on \J l^'u) A for various Ra. a) semi-log plot and non- 
dimensionalized with friction velocity Ut. Dashed line rep- 
resents — z'^; b) normalized by the maximum velocity 
timax and the kinematic boundary layer thickness A„. 



Shown in Fig. ^Bp is the velocity profile normalized 
by the outer variables, i.e. the boundary layer thickness 
An and the maximum velocity Wmax = u{\u)- Although 
the profiles show that there is a Ra dependence, it is 
very weak. The weak influence of the Ra number - es- 
pecially for the two lower Ra numbers considered - is 
further evidence that the kinematic boundary layer does 
not behave as a classical forced turbulent boundary layer. 
Note that the approximate universality of the velocity 
profiles means that inner and outer variables can be in- 
terchanged, in the sense that dzu\yj oc Umax/Au. 

Several experiments have shown universality in Ra 
upon an outer scaling by boundary layer thickness and 
maximum velocity [J [la . , so it is quite interesting 
that the boundary layer profile found here (Fig.[2}3) has a 
(small) Ra dependence. There may be several reasons for 
this difference. The experiments have been carried out 
at much higher Ra, in the range Ra = 2 x 10® — 9 x 10^ 
and at higher Pr (the working fluid was water). Further- 
more, the presence of side walls and the smaller aspect 
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TABLE I: Characteristic numbers for the boundary layer 
profile at various Ra: the shear- Reynolds number Re^, the 
shape factor S and the friction coefficient C/. 



Ra 


Re^ 


Cf 


S 


1.15 X fO^ 


26 


1.02 


2.37 


f .00 X fO'' 


52 


0.51 


2.35 


1.00 X 10^ 


119 


0.23 


2.27 



ratio will be of influence. 

It is useful to express the shear-Reynolds number in 
terms of Re, A„ and the non-dimensional velocity gradi- 
ent at the wall. Let the outer scaled variables be denoted 
by z = z/Xu and u = u/umax- The non-dimensional ve- 
locity gradient at the wall is connected to the wall-shear 
stress by r^, = pvUraaxX~^ dzu\'un where dzu\w is the 
non-dimensional velocity gradient at the wall. Hence, 
the shear Reynolds number can be expressed as 



Rc^ = Re 



(6) 



All three terms Re, and dzu\w depend on Ra, al- 
though the Ra dependence of the last term is very weak 
as dzu\w (X Ra^"®. 



C. Friction and shape factor 

The friction and shape factor [e.g. @, have been 
calculated for all three Ra (Table [J). The friction factor 
is defined as 



C 



.Ref 



Rp2 
2 "max '-^'^ 



(7) 



Here we note that combining ([6]), ([7|) and neglecting the 
small Ra dependence of the wall-gradient dzu\w gives 
that Cf can be approximated by 
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FIG. 3: y-averaged wind structure at Ra — 10® and Pr = 1. 
Colorscheme is by the relative temperature Qr = Q — {Q)a- 
Dark areas are relatively cold, and white areas are relatively 
hot. 



plate flow Cf ~ 0.664Re^ ' . Hence, judging from the 
scaling of friction factor, the behavior of the boundary 
layer would be classified as laminar. These results are 
consistent with 0, El, El • 

The shape-factor S is defined as 5 = Si/ 62, where i5i 
and S2 are the displacement and momentum thickness, 
given by: 



Si^ (1 - 



-(1 



-)dz, 



-)dz. 



S2 



For laminar profiles, such as Poiseuille flow and the Bla- 
sius solution for the developing flow over a flat plate, 
the shape factor is approximately 2.5 [e.g. @, [l^. For 
turbulent plane channel flow, flat-plate constant-pressure 
boundary layers and a plane turbulent wall jet jT9| the 
shape factor is approximately 1.3 — 1.4. Based on this 
information, the values from Table U indicate that the 
velocity profile follows a laminar-like distribution with a 
slight trend towards turbulent values as Ra increases. 

If the shape and friction factor are taken to be rep- 
resentative to distinguish a laminar from a turbulent 
boundary layer, the boundary layer would be classified 
as laminar. In the next sections we will study the mo- 
mentum budgets of the boundary layers, and compare the 
time and space spectra of boundary layer and the bulk. 
It will be shown that from this perspective, the kinematic 
boundary layer has many features of turbulence. 



r ^ 2 fXu 



(8) 



This is consistent with the approximation r^, « vU/Xu, 
which is an important assumption in the Grossmann- 
Lohse theory [l||. The observation ^ will prove to be 
important to establish the scaling of Xu/H in section HVl 
Based on the values of Table U and in terms of Re, 
the friction factor Cf scales as C/ cx Re^'''®*'. An 
empirical relation for turbulent plane channel flow is 
Cf = 0.073Re~°'^'^, with Re based on channel half width 
and mean velocity across the channel [l^. The friction 
factor of laminar boundary layers have a stronger de- 
pendence on Re; for plane Poiseuille flow Cf = 8/Re 
(Re based on full channel height) and for the Blasius flat 



D. Fluctuations and spectra 



In Fig. [3^-c the average velocity proflle (M)j//'u,„ax is 
shown for Ra = 1.15 x 10^, 10^ and 10'', together with the 
turbulence intensity of the horizontal and vertical fluctu- 



atlons, (u'u 



/"max respectively. 
These are the profiles of the j/-avcraged wind structure 
(see Fig. [3]), with the x-location chosen such that the 
horizontal velocity is at its maximum, i.e. where the flow 
is parallel to the wall and from left to right. A strik- 
ing feature of the turbulence intensity of the horizontal 
fluctuations, is that it is so large compared to the mean 
wind, namely 70-80%. For turbulent channel flow, typi- 
cal turbulence intensities are 5-10%. Outside the thermal 
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FIG. 4: Close-up of horizontal velocity and turbulence inten- 
sities (legend in Fig. |4Jd) at the position with the maximum 
horizontal velocity, (a) Ra = 1.15 x lO''; (b) Ra = 1.00 x 10*^; 
(c) Ra = 1.00 X 10^. The horizontal dashed line indicates the 
edge of the thermal boundary layer Ae. 



boundary layer the horizontal turbulence intensity is con- 
stant. The vertical turbulence intensity is not as large as 
the horizontal due to wall blocking, but is still 20% at the 
edge of the thermal boundary layer, and 50% at the edge 
of the kinematic boundary layer. This confirms that fluc- 
tuations in large aspect-ratio domains are larger relative 
to the wind [20|. in comparison with small aspect-ratio 
domains (e.g. reports turbulence intensities of 20%). 

One of the main features of turbulence is the presence 
of a continuous range of active scales. A simulation at 
Ra = 10^ is used to obtain both spatial and temporal 
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FIG. 5: Spectra of the horizontal velocity components in the 
boundary layer and in the bulk at Ra = 10^ normalized by the 
bulk variance, (a) temporal spectrum (b) spatial spectrum. 



spectra of the horizontal velocity components. To collect 
temporal spectra, eight points have been monitored: four 
bulk and four boundary layer points. The bulk points are 
taken at Zbuik = Hj^ and the boundary layer points were 
chosen according to Zbi = A„. The temporal spectra are 
generated by segmenting the time series and a Welch win- 
dow has been used. Then, averaging was performed over 
the spectra of the two horizontal velocity components 
and the four monitoring points. The spatial spectra were 
collected by performing a 2D FFT and integrating over 
circles k'^ + k^ = k^ and averaging over approximately 10 
turnovers. 

The temporal spectra of the horizontal velocity com- 
ponents at Ra = 10^ are shown in Fig. [5^. There is a 
continuous range of active scales which spans about two 
decades, although turbulence production and dissipation 
are not sufRciently separated to form a clearly discernible 
inertial subrange. The spatial spectra (Fig. \5]p) also re- 
veal a continuous range of active scales. 

What is striking about the spectra of the bulk and the 
boundary layer is how similar they are, both in range of 
active scales and in amplitude. Despite a mild damping 
at the intermediate frequencies and wave numbers, the 
similarity indicates that the dynamics of the bulk and 
the boundary layer - both temporal and spatial - are 
very similar. We note that the simulation at Ra = 10^ 
is well inside the hard-turbulence regime. The transition 
to hard turbulence occurs at much lower Ra for large 
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TABLE II: Budget terms for momentum and heat equation. 
A V P B 7^ 

dtUi = —djUjUi +ud^Ui —dip +/3g&5is —dju'^u'- 
dtQ= -djUjQ +K9|e -dj^' 



aspect-ratio domains than the generally accepted value 
of Ra = 4 X 10^ [llj. Indeed, for aspect-ratio 6 domains, 
the flow has hard-turbulence regime scaling occurs from 
Ra = 5 X 10** upwards [1]. If one accepts that the flow 
core is turbulent, then Fig.[5]indicates that the boundary 
layers are turbulent as well. 

The striking similarity between the spectra in the bulk 
and the boundary layers seems to be a robust and general 
feature of Rayleigh-Benard convection. In a recent paper 
[l3l |. we present combined experimental and numerical 
results of an aspect-ratio 4 cavity filled with water for 
Rayleigh numbers ranging from 5 x 10'' to 10^. For all 
Ra from 10^ upwards, it is found that the spectra in the 
bulk and the boundary layer are practically identical. 



E. Momentum budgets 

Momentum budgets are a very direct way to get an 
impression of the importance of the turbulent Reynolds 
stresses. As before, y-averaged results (Fig. [3]) are used 
for convenience of presentation. Checks have been made 
to ensure that the budgets shown here are also repre- 
sentative for the three-dimensional field. The a;-location 
has been chosen such that the horizontal velocity is at 
its maximum, i.e. where the flow is parallel to the wall 
and from left to right. This guarantees that horizon- 
tal gradients are small, and that no adverse of favorable 
pressure gradients are present. Shown are budgets for 
Ra = 1.15 x 10^ (Fig[nii,d,g), Ra = 1.00 x 10^ (Fig[ni3,e,h) 
and Ra = 1.00 x lO"^ (Fig[6t,f,i). The budgets for the hor- 
izontal (Figs, in^-c) and vertical momentum (Figs. [SJl-f) 
have been nondimensionalized by / H ~ PgAQ, and 
heat budget ^-i) by AeU/H ^ ^l3g{AQYlH. The 
legend for the budgets is shown in Fig. [6^; and the bud- 
get terms are defined in Table Ull The z-coordinate has 
been scaled by A„ and the horizontal dashed line denotes 
z = Ae- For reference, the ratio Ae/A„ is 0.8, 0.6 and 
0.38 for the simulations at Ra = 10^, 10^ and 10^ respec- 
tively. 

For the horizontal momentum budgets (Figs.[6^-c), the 
balance is between the horizontal pressure gradient V 
and diffusion T) for z < Ae ■ Outside the thermal bound- 
ary layer, TZ is not negligible; on the contrary, TZ fully 
balances the pressure gradient V near z = A„. This indi- 
cates that the turbulence outside the thermal boundary 
layer is key to the boundary layer thickness, as will be 
outlined in section IIVI As the location of the budgets 



has been chosen such that all horizontal derivatives are 
small, V K, vdlu and TZ k, —dzw'u'. 

Log-scaling is expected in the inner layer where w'u' is 
constant, so that TZ = —dzW'u' = 0. For channel flow, TZ 
is zero at the wall and peaks in the buffer layer, marking 
the transport of momentum from the outer to the inner 
layer. After the peak, it crosses the zero axis where the 
log-layer is expected. This behavior of TZ is absent for 
Ra = 1.15 X 10^, but as Ra increases a peak forms in- 
side the thermal boundary layer (Fig.[6K-c). However, in 
terms of forcing V is always much larger than the small 
peak TZ for the range of Ra under consideration, which 
again confirms that this is not a classical forced turbulent 
boundary layer. 

Figs. [6jl-f show the budgets of the w-momentum equa- 
tion. Here the balance is between buoyancy B, the 
vertical pressure gradient T^ and the Reynolds stresses 
TZ ~ —dzw'w'. Very near the wall, roughly in the lower 
half of the thermal boundary layer, the buoyancy B and 
pressure T' are in balance, so the flow is neutrally buoy- 
ant here. Further away from the wall, at the edge of the 
thermal boundary layer, the contribution of TZ is signifi- 
cant, even if it may seem small compared to the near- wall 
(hydrostatic) balance of T^ and B. In fact, comparing TZ 
of the vertical momentum equation to the magnitude of 
terms in the horizontal momentum equation shows that it 
is of the same magnitude as —dxP- Outside the boundary 
layer, the pressure gradient P is positive and is balanced 
purely by fluctuations TZ. 

The ©-momentum budgets (Figs. [Gji-i) show a bal- 
ance between thermal diffusion 2? = kO^Q, turbulence 
TZ = —dzw'Q' and there is a contribution from advcction 
A. Judging from the peak of A around z/A^ = 1, the 
nonzero contribution of A to the heat budget is prob- 
ably caused by some spatial variations in O by which 
dxuQ 0. The peak of T> and TZ is always located just 
inside the thermal boundary layer, representing the lo- 
cation where diffusion and fluctuations most effectively 
exchange heat. 

It is striking that the dominant length scale for the 
budgets is the thermal boundary layer thickness Ae 
(which is denoted by the horizontal dashed line in Fig.lH]), 
and not as one may expect, the kinematic boundary layer 
thickness. Perhaps this should not be too much of a sur- 
prise, as the thermal boundary layer thickness can be well 
represented by Ae = -ff/(2Nu), and the Nusselt number 
Nu represents the efficiency of the convective heat trans- 
fer mechanism of the flow, resulting from the non-linear 
coupling of temperature and velocity under the action 
of buoyancy. Therefore, Ae is equally important for the 
heat-budget and for the momentum budgets. In fact, Ae 
is a dominant parameter in the scaling of both C/ and 
A,j, as will be shown in section Hvl 

The findings of Figs. [S^-i can be summarized as follows 
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FIG. 6: Momentum- and heat-budget in the boundary layer at the position with the maximum horizontal velocity for various 
Ra: (a,d,g): Ra = 1.15 x 10^ (b,e,h): Ra = 1.00 x 10*^; (c,f,i): Ra = 1.00 x 10^.. (a,b,c); S-momentum budget; (d,e,f): 
w-momentum budget; (g,h,i): 0-budget. The legend for Figs, (a-i) is shown in Fig. [6J5 and the horizontal dashed line indicates 
the edge of the thermal boundary layer z = \q. 



for the u, w and budget, respectively: 

dxP + dzw'u' ~ vd^u, 
dzP + dzw'w' ~ PgQ, 



(9) 
(10) 
(11) 



These equations represent the boundary layer equations 
at the x-location where the flow is parallel to the wall and 



horizontal derivatives are negligible (roughly halfway be- 
tween the impingement and detachment region). Note 
that even though the ui-momentum equation is not di- 
rectly coupled to the other two equations, the vertical 
fluctuations w'w' are non-trivially coupled to w'u' and 
w'Q' as these terms represent to a large extent the plumes 
emerging from and impinging on the boundary layers. 
The equations above are two-dimensional, but by the ab- 
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TABLE III: Decomposition of the friction factor C/ according 
to JH]) 



Ra 


Cf = 






10^ 


0.81 


0.88 


-0.07 


10" 


0.39 


0.42 


-0.03 


10^ 


0.17 


0.19 


-0.02 



sence of transversal derivatives, it can be expected that 
these equations are valid for the three-dimensional case 
as well, in a local coordinate system aligned with the flow 
and at the location where the flow is parallel to the wall. 

The boundary layer equation ^ clearly shows that one 
cannot neglect the influence of turbulence in the bound- 
ary layer dynamics. Hence, the laminar boundary layer 
equation ([3]), which lead to the scaling Au oc Re~^^^ is 
not valid: additional information is required about w'u' 
to estimate A^. In section ITVl the scaling behavior of A„ 
will be derived using flow-specific information obtained 
from the DNS results. 



F. The friction factor decomposed 

By using the boundary layer equation Q, the domi- 
nant contributor to the friction factor can be identified. 
Integrating Q over the kinematic boundary layer and 
substituting ([7]), the friction factor Cf is composed of a 
contribution from pressure and a turbulent momentum 
flux as 



9l 

2 



1 



{V + n)dz 



dxpdz 



w'u'\\^ 



(12) 



The terms on the right hand side of have been cal- 
culated with the DNS results and are presented in Table 
|TIT][25| The decomposition clearly demonstrates that Cf 
is dominated by the pressure gradient. The turbulent 
momentum flux w'u' is small but positive, i.e. a flux out 
of the boundary layer. Hence, we conclude that the dy- 
namics of the wall-friction is not governed by turbulence 
(i.e. Reynolds stress) as in a forced turbulent boundary 
layer. In the latter, the dynamics is dominated by a large 
momentum-flux into the boundary layer, while the con- 
tribution of the pressure gradient is negligible. 

Using p^ . Cf can be parameterised. Shown in Fig. 
[7] is the effective forcing 7^ + 71 = —dxp — dzw'u'. To 
first order, for z < Xq, V + TZ ~ V while for z > Ae 
the Reynolds stress forcing TZ balances the pressure V to 
that V + TZ 0. Hence, Cf can be estimated by 




0.025 0.05 0.075 
V + TZ 



0.1 



FIG. 7: The effective forcing 1^ + 71= —dxP — dzW'u', which 
vanishes quickly outside the thermal boundary layer. 



Cf 



Ae 



Vdz : 



2Ae 



dxp\. 



(13) 



Clearly, (fT3|) holds at moderate Ra only, when turbu- 
lent shear production in the boundary layer is small. The 
formation of the peak inside the thermal boundary layer 
at Ra = 10^ and Ra = 10'' (Fig. [7|) suggests that shear 
production becomes more important as Ra increases, and 
this will have to be accounted for in ([T3|) at higher Ra. 
However, it was shown in the accompanying paper Q 
that the wind velocity becomes independent of Cf at suf- 
ficiently high Ra, because Cf is negligible compared to 
the mixing parameter a. Therefore, incorrect scaling be- 
havior in (|13p will not influence the wind dynamics at 
high Ra. 



IV. SCALING OF Cf AND A„ 

Using the simple two-equation wind model derived in 
the accompanying paper [9| , we can establish the scaling 
behavior of C/ and A^, . The model uses a dimension- 
less wind- velocity Uw ~ U^/Uf and spatial temperature 
difference Qw = Q^/AQ, where Uf = y/JJgAQH is the 
free-fall velocity. The governing equations of the model 
are given by 



2Ll 



1 



dt 2Ll + 1 \2L. 
dQi^ 2Ae 4q! 



■e^-(4a + C/) 



U,, 



(14) 



dt 



LlPn 



Uu 



AeRe/Pr 



(15) 



Here ~ Ly,/ H where Lw is the typical roll size, Ae = 
Xe/H, A„ = Xu/H and Re/ = UfH/f = Ra^/^Pr"^/^ 
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The turbulent Prandtl number Pr^ and the mixing pa- 
rameter a are coefficients with values 0.85 @ and 0.6 re- 
spectively. The pressure difference which drives the wind 
is generated by a spatial temperature difference Qw (it is 
relatively hot where the flow ascends and relatively cold 
where it descends, see Fig. [3]). The temperature differ- 
ence is in its turn generated by large horizontal heat 
fluxes originating from the interaction between the mean 
wind and temperature field. The model depends on Ra, 
Pr and i^,, where Ae = Ae(Ra, Pr) and C/ = C/(Ra, Pr) 
have to be provided. Based on the analysis of the friction 
factor fsection flllFj) . an explicit expression for Cf can be 
derived, by which the model only depends on empirical 
input for Ae (and thus Nu). 

The steady state estimate for the pressure gradient at 
the bottom wall of the wind model is Q 



dxP\, 



PgH 



Using (fT6|) . Cf (fT3|) can be further specified as 



C 



2Ac- 



/ 



2Ae e 



H Lu, Ul Ae 



LwU^ 



(16) 



(17) 



Hence, the wall friction term is linear in the temperature 
difference as 



Cf 



2Ar- 



(18) 



Here, we assumed that sgn Uw = sgn Q^. With p8|) . the 
empirical specification of C/(Ra, Pr) is no longer neces- 
sary, and the model is given by: 



dUw 
di 

di 



2L, 



2LI 
2Ae 



1 - 4Ae ^ 



u„ 



2L^ 
4a 



4a 













2 


Ae 


Rc/P 



(19) 



(20) 



The steady state solution of the model as a function 
of Ra is shown in Fig. [51 At this point, the only empir- 
ical data used in the model is the powerlaw for Ae and 
the roll size L^. The mixing parameter a is kept at the 
same value as in [i], namely 0.6. As can be seen, the 
model captures the trends of U^,, ©tu and Cf satisfac- 
torily. Note that the profiles could be made to match 
quantitatively as well when one would introduce some 
additional coefficients. However, the focus of this paper 
is not to develop a carefully tuned model, but to elicit 
general scaling behavior. 

It is not very useful to have an expression for Cf in 
terms of ©„,, as this quantity is rarely measured. How- 
ever, by using the steady state solution of p9)) . G) can be 
expressed in terms of Uw as 



SaLu, 
1-4A0 



(21) 



so that the 8^: dependence of Cf can be eliminated. Us- 
ing HH]) and C/ is given by 



Cf 



16a 



1 - 4Af- 



-Ar- 



(22) 



Hence, when Ae 1, the model predicts that C/ oc Ae- 
Note that the Ae term in the denominator represents the 
effects of wall friction. Hence, when Ae ^ 1, C f scales 
independently of wall-effects. It is the turbulence in the 
outer flow which fully determines the velocity at the edge 
of the boundary layer. 

A scaling relation for A„ can be derived by using the 
two different expressions for C/, ([8]) and (fT7|) . This re- 
sults in 



2Re" 



2Ar- 



©„ 



Using Re 



Re/, A„ is approximated by 



All 



Lw 






Ae 


On, 


Re/ 



(23) 



Dropping the absolute signs and using (j2ip . A„ is given 

by 



Au 



l-4Af 



4Ae 1 



SaAeC^toRe/ 



8a AeRe 



(24) 



Upon assuming that Ae <C 1 , it follows that A„ scales as 



Au (X Aq Re" 



^ . Fig. [8jl shows the prediction of the wind 
Although the boundary layer thickness is 



model for A 

underprcdicted, the trend is in agreement with the DNS 
data. Given the simplicity (with only one calibration 
parameter a), the model captures the trends of wind ve- 
locity, spatial temperature difference, friction factor and 
kinematic boundary layer thickness satisfactorily. 



V. TURBULENT OR NOT? 

The apparently contradicting findings reported in the 
previous sections is quite intriguing. On the one hand, 
the results indicate that the kinematic boundary layer 
is turbulent. The deduced boundary layer equation (O 
shows that forcing due to turbulent Reynolds stresses 
is significant, in particular outside the thermal bound- 
ary layer. Furthermore, the spectra in the bulk and the 
boundary layers are nearly indistinguishable and show 
the existence of a continuous range of active scales both 
in space and time. Both are an indication for turbulence. 

On the other hand, the results suggest that the kine- 
matic boundary layer does not correspond to a classi- 
cal turbulent boundary layer. The Reynolds numbers 
in the Ra range we consider (Re ~ 1500 at Ra = 10^) 
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FIG. 8: The predictions of the wind model (eqns p9p . (|20|l . thick lines) compared to the DNS results (diamonds) for a) the 
typical wind Uw', b) the spatial temperature difference Qw', c) the friction factor C/ and d) kinematic boundary layer thickness 



are generally considered too low to sustain turbulence. 
Moreover, the friction factor C/ for a classical forced 
boundary layer has a weak dependence on Re (reflect- 
ing the quadratic wall friction), and is dominated by the 
turbulent momentum-flux from the free stream. For the 
boundary layers under consideration, the dominant con- 
tributor to C/ is the pressure gradient fsection rill F[) and 
not the momentum- flux. Consequently, C/ has a signifi- 
cant Re dependence. The near-universal profiles (found 
in the present work especially for the two lower Ra num- 
bers) as a function of Ra based on the outer variables A„ 
and Umax (Fig- [l]) arc further evidence against a classi- 
cal turbulent boundary layer: a turbulent boundary layer 
can by definition not be universally scaled by outer vari- 
ables. 

The difference between classical forced turbulence 
boundary layers and a boundary layer of Raylcigh- 
Benard convection may be best characterized by the way 
the turbulence is produced and redistributed. For forced 
flow cases, turbulence cannot be maintained at low Re, 
as the dissipation in the boundary layer will be stronger 
than the production. However, for Rayleigh-Benard con- 
vection the production and transport of turbulent kinetic 
energy (TKE) is not confined to the inner layer alone. In- 
stead, TKE is produced in the bulk, where it is partially 



dissipated. The surplus is transported to the bound- 
ary layer by pressure velocity fluctuations [see also [2^ . 
Therefore, there is no local equilibrium between produc- 
tion and dissipation, and turbulence can be maintained 
in the boundary layers below the critical Re. At suf- 
ficiently high Ra, instabilities due to shear can be ex- 
pected to maintain themselves, and several experiments 
and simulations show such a transition around Ra = 10^^ 

A simple explanation for the laminar-like scaling of 
classical integral boundary layer parameters may be that 
the forcing in the wall-parallel direction is very weak com- 
pared to the forcing in the wall-normal direction (plume 
impingement and detachment). Indeed, the forcing in the 
vertical direction is the direct result of buoyancy, while 
the pressure gradient in the horizontal direction forms is 
due to large-scale differences in mean temperature. This 
can be made explicit by considering the ratio of forcing 
in the wall-normal direction (buoyancy) and wall-parallel 
direction (|16|) . which is given by 



/?5Ae 



= Lu.e: 



(25) 



At Ra =10^, this ratio is approximately 50, and at Ra = 
10^, the ratio is approximately 100. Thus, the boundary 
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layers under consideration here are forced primarily in 
the wall-normal direction, and the force generating the 
wind is relatively weak. 

Despite the laminar-like scaling of the integral parame- 
ters, a parallel can be drawn with a fully developed forced 
boundary layer: both have a viscous sublayer dominated 
by viscosity which suppresses instabilities and prevents 
their growth and development of turbulence. However, 
as demonstrated by seminal experiments in the sixties 
[23 |. despite linear velocity variation, the flow within the 
sublayer in a forced boundary layer is not laminar, but ac- 
companied by considerable irregular fluctuations, streaks 
and other structures. One can argue that the same dy- 
namics occur in the boundary layers of Rayleigh-Benard 
convection. In particular. Figs. [6K-c indicate that the 
thermal boundary layer < z < Ae functions as a vis- 
cous sublayer, and the region Ae < ^ < A„ as an cross- 
over region between the exterior flow and the thermal 
boundary layer. The absence of a constant stress layer 
dominated by the turbulent momentum flux suppresses 
a logarithmic region and marks a fundamental difference 
with forced turbulent boundary layers. 

VI. CONCLUSIONS 

The aim of this paper has been to study the bound- 
ary layers which develop under the joint action of plumes 
and wind in Rayleigh-Benard convection at the top and 
bottom plates. Direct numerical simulation was used for 
simulations at Ra = {10^ 10^ 10^ 10^} and Pr = 1 for 
r — 4 aspect ratio domains with periodic side bound- 
ary conditions. For each Ra, 10 independent simulations 
have been carried out, resulting in approximately 400 
independent realizations per Ra. Processing the results 
using symmetry-accounting ensemble averaging made it 
possible to retain the wind structure, which would nor- 
mally cancel out due to the translational invariance of 
the system. 

The importance of Reynolds-stresses in the boundary 



layers, as well as the temporal and spatial spectra in- 
dicate undoubtedly a turbulent character of the bound- 
ary layer. However, the behavior is rather different from 
classical forced boundary layers, as can be judged from 
the laminar-like scaling of the classical integral boundary 
layer parameters. Indeed, viscous effects play an impor- 
tant role within the thermal boundary layer, and a large 
turbulent momentum-flux from the external stream is ab- 
sent. This difference is probably caused by the fact that 
the turbulence inside the kinematic boundary layer of RB 
originates from the bulk, whereas classical forced bound- 
ary layers are in a local equilibrium between production 
and dissipation of turbulent kinetic energy. 

Due to the importance of Reynolds stresses in the 
boundary layer, the arguments underpinning the kine- 
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matic boundary layer scaling A„ cx Re~"^^ do not hold. 
Using the DNS results and a conceptual wind model 0], 
explicit expressions for C/ and A„ were derived. It was 
found that the friction factor should scale proportional to 
the thermal boundary layer thickness as C/ cx Ae. The 
kinematic boundary layer thickness A„ scales inversely 
proportional to the thermal boundary layer thickness and 
the Reynolds number as A„ cx Re~^AQ^. The predicted 
trends for C/ and Au are in agreement with the DNS 
results. 

With the closure for C/, the model (^0]) depends 
solely on empirical input for Ae, and predicts the wind 
Reynolds number Re, friction factor Cj and kinematic 
boundary layer thickness A„. 



Acknowledgments 

This work is part of the research programme of the 
Stichting voor Fundamenteel Onderzoek der Materie 
(FOM), which is financially supported by the Ned- 
crlandsc Organisatic voor Wctenschappclijk Onderzoek 
(NWO). The computations were sponsored by the Sticht- 
ing Nationalc Computerfacilitcitcn (NCF). 



[1] S. Grossmann and D. Lohse, J. Fluid Mech. 407, 27 
(2000). 

[2] Y. B. Xin, K. Q. Xia, and P. Tong, Phys. Rev. Lett. 77, 
1266 (1996). 

[3] Y. B. Xin and K. Q. Xia, Phys. Rev. E 56, 3010 (1997). 
[4] R. M. Kerr, J. Fluid Mech. 310, 139 (1996). 
[5] R. M. Kerr and J. R. Herring, J. Fluid Mech. 419, 325 
(2000). 

[6] H. Schlichting and K. Gersten, Boundary layer theory 

(McGraw-Hill, 2000). 
[7] S. Grossmann and D. Lohse, Phys. Fluids 16, 4462 

(2004). 

[8] S. Grossmann and D. Lohse, J. Fluid Mech. 486, 105 
(2003). 

[9] M. van Reeuwijk, H. J. J. Jonker, and K. Han- 
jalic, Submitted to Phys. Rev. E (2007), URL 



'http : //arxiv . org/abs/0709 . 0304 ' 

[10] X. Chavannc, F. Chilla, B. Castaing, B. Hebral, 
B. Chabaud, and J. Chaussy, Phys. Rev. Lett. 79, 3648 
(1997). 

[11] X. Chavanne, F. Chilla, B. Chabaud, B. Castaing, and 

B. Hebral, Phys. Fluids 13, 1300 (2001). 
[12] G. Amati, K. Koal, F. Massaioli, K. R. Sreenivasan, and 

R. Verzicco, Phys. Fluids 17, 121701 (2005). 
[13] J. Verdoold, M. van Reeuwijk, M. J. Tummers, H. J. J. 

Jonker, and K. Hanjahc, Phys. Rev. E 77, 016303 (2008), 

URL http : //ar xiv . org/abs/0707 . 2485 , 
[14] M. van Reeuwijk, H. J. J. Jonker, and K. Hanjalic, Phys. 

Fluids 17, 051704 (2005). 
[15] F. M. White, Viscous fluid flow (McGraw-Hill, 1991). 
[16] S. Lam, X. D. Shang, S. Q. Zhou, and K. Q. Xia, Phys. 

Rev. E 65, 066306 (2002). 



12 



[17] X. L. Qiu and K. Q. Xia, Phys. Rev. E 58, 5816 (1998). 
[18] R. B. Dean, J. Fluid Eng.-T. ASME 100, 215 (1978). 
[19] N. Rajaratnam, Turbulent jets, no. 5 in Developments in 

water science (Elsevier, 1976). 
[20] J. J. Niemela and K. R. Sreenivasan, J. Fluid Mech. 557, 

411 (2006). 

[21] F. Heslot, B. Castaing, and A. Libchaber, Phys. Rev. A 

36, 5870 (1987). 
[22] R. M. Kerr, Phys. Rev. Lett. 87, 244502 (2001). 



[23] J. J. Niemela and K. R. Sreenivasan, J. Fluid Mech. 481, 
355 (2003). 

[24] S. J. Kline, W. C. Reynolds, F. A. Schraub, and P. W. 
Runstadler, J. Fluid Mech. 30, 741 (1967). 

[25] The values for C/ in Table IIIII are slightly smaller than 
those in Table U] as the latter is deduced from the plane- 
averaged squared mean velocity. 



